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[i]  A  numerical  simulation  of  the  low-altitude  electron  radiation  belt  is  described.  It 
includes  dependences  on  the  electrons’  bounce  and  drift  phases,  equatorial  pitch  angle, 
and  kinetic  energy  in  the  range  of  ~1  to  several  MeV  at  L  =  3.5.  Physical  processes  in 
addition  to  the  adiabatic  electron  motion  are  pitch  angle  diffusion  and  backscattering  from 
a  realistic  model  atmosphere.  Quasi-linear  diffusion  coefficients  are  calculated  from  a 
model  of  the  whistler  mode  plasmaspheric  hiss  wave  intensity.  Comparisons  of  the 
simulation  results  with  electron  data  from  a  low-altitude  satellite  show  that  the  model 
accounts  for  the  main  features  of  the  electron  spatial  distribution  during  selected  periods  of 
differing  geomagnetic  activity.  INDEX  TERMS:  2720  Magnetospheric  Physics:  Energetic  particles, 
trapped;  2716  Magnetospheric  Physics:  Energetic  particles,  precipitating;  2753  Magnetospheric  Physics; 

Numerical  modeling;  2772  Magnetospheric  Physics:  Plasma  waves  and  instabilities;  KEYWORDS:  radiation 
belt,  electrons,  diffusion,  atmospheric  scattering 

Citation:  Selesnick,  R.  S.,  M.  D.  Looper,  and  J.  M.  Albert  (2004),  Low-altitude  distribution  of  radiation  belt  electrons,  J.  Geophys. 
Res.,  109,  A11209,  doi:10.1029/2004JA01061 1. 


1.  Introduction 

[2]  The  distribution  of  radiation  belt  electrons  is  consid¬ 
erably  more  complex  at  low  altitudes  than  at  high  altitudes 
because  of  the  interaction  with  the  atmosphere  combined 
with  the  non-dipolar  nature  of  the  geomagnetic  field. 
Radiation  belt  populations  whose  mirror  points  do  not  reach 
the  atmosphere  can  usually  be  described  by  an  average  over 
all  three  phases  of  the  adiabatic  particle  motion;  the  gyration 
phase  about  the  local  magnetic  field  direction,  the  bounce 
phase  along  the  magnetic  field  line  between  the  mirror 
points,  and  the  drift  phase  around  the  Earth.  This  averaging 
reduces  the  dimensionality  of  the  phase  space  density  /from 
6  to  3.  The  asymmetric  geomagnetic  field  allows  the 
possibility  of  mirror  points  both  above  and  below  the  top 
of  the  atmosphere  within  the  adiabatic  motion  of  a  single 
particle.  In  a  population  where  a  substantial  fraction  of  the 
mirror  points  are  within  or  below  the  dense  atmosphere, /is 
likely  to  have  a  strong  dependence  on  the  both  the  bounce 
and  drift  phases,  so  that  averaging  over  those  two  variables  is 
inadvisable  and  the  dimensionality  of /is  reduced  then  only 
to  5.  Satellites  in  low  Earth  orbit  (with  altitudes  103  km) 
encounter  such  a  population  of  radiation  belt  electrons. 
Therefore  it  is  desirable  that  radiation  belt  models  for 
application  to  such  orbits  should  retain  both  the  bounce 
and  drift  phase  dependencies. 

[3]  Low-altitude  satellites  spend  a  large  fraction  of  their 
observing  time  within  the  bounce  loss  cone,  that  region  of 
momentum  space  for  which  the  mirror  points  are  below  the 
atmosphere  in  either  hemisphere.  Numerical  modeling  of 
the  bounce  loss  cone  electron  distribution,  including  the 
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bounce  phase  dependence  and  electron  backscattering  from 
the  atmosphere,  was  described  by  Davidson  and  Walt 
[1977],  An  analytic  approximation  with  asymptotic  match¬ 
ing  of  the  trapped  electron  distribution  was  described  by 
Retterer  et  al.  [1983].  These  models  were  applicable  to  the 
auroral  electron  kinetic  energies  of  ~10  keV,  for  which  the 
drift  timescales  are  relatively  long  and  the  drift  phase 
dependence  of  the  electron  distribution  can  reasonably  be 
neglected. 

[4]  For  electrons  with  kinetic  energies  ~1  MeV  the 
azimuthal  drift  motion  is  significant.  A  numerical  model  of 
their  distribution,  including  the  drift  phase  dependence  and 
pitch  angle  diffusion  from  the  stably  trapped  population  into 
the  quasi-trapped  population  was  described  by  Selesnick  et 
al.  [2003],  As  the  geographic  longitude  of  a  low-altitude 
satellite  varies,  either  by  its  orbital  motion  or  by  the  rotation 
of  the  Earth,  it  is  accessible  to  a  varying  fraction  of  electron 
trajectories  that  are  either  stably  trapped,  with  the  drift  shell 
entirely  above  the  dense  atmosphere,  or  quasi-trapped,  with 
mirror  points  below  the  dense  atmosphere  at  a  remote 
longitude.  The  quasi-trapped  electrons  are  outside  the  local 
bounce  loss  cone  but  are  in  the  drift  loss  cone  because  their 
drift  motion  will  eventually  take  them  into  the  bounce  loss 
cone.  The  trapped  and  quasi-trapped  electron  intensities  are 
usually  significantly  higher  than  those  of  the  untrapped  (or 
bounce  loss  cone)  electrons,  so  they  can  be  the  dominant 
contributors  to  the  total  electron  fluence  over  an  orbital 
timescale. 

[5]  The  energetic  ~1  MeV  electrons  can  also  include  a 
significant  bounce  loss  cone  population.  To  describe  this 
population  it  is  necessary  to  combine  the  two  types  of 
calculations  described  above  into  a  single  model  that 
includes  both  the  bounce  and  drift  phase  dependences.  That 
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is  the  goal  of  this  work.  The  model  has  three  main  inputs:  a 
description  of  the  electron  scattering  by  the  atmosphere 
based  on  Monte  Carlo  simulations,  the  pitch  angle  diffusion 
coefficients  derived  from  the  theory  of  wave  particle  inter¬ 
actions,  and  the  initial  stably  trapped  electron  distribution 
obtained  from  satellite  measurements.  By  comparing  the 
model  results  with  observations  of  the  bounce  and  drift  loss 
cone  electron  distributions  we  can  then  determine  whether 
these  model  elements  are  sufficient  to  describe  the  main 
spatial  variations  in  those  electron  populations.  Such  calcu¬ 
lations  can  also  be  of  practical  value  for  predicting  the  range 
of  electron  intensities  in  various  low  Earth  orbits. 


2.  Model  Description 
2.1.  Transport  Equation 

[6]  Stochastic  pitch  angle  scattering  at  constant  energy 
can  be  described  by  a  one-dimensional  Fokker-Plank  equa¬ 
tion  based  on  the  spherical  volume  element  in  momentum 
space,  with  only  the  pitch  angle  a  allowed  to  vary.  With  the 
condition  that  the  steady  state  solution  be  isotropic,  the 
friction  and  diffusion  terms  can  be  combined.  The  resulting 
diffusion  equation  is 


df 

dt 


1  d 
sin  a.  dot 


fsinaZ>aci 


(1) 


where  6(5)  is  the  Dirac  delta  function.  All  electrons  reaching 
the  location  of  the  dense  atmosphere,  at  s  =  s%m  or  in 
the  northern  or  southern  hemisphere  respectively,  are 
considered  to  be  lost.  The  corresponding  backscattered 
distributions  at  the  same  locations,  fis  and  f£s,  which 
depend  on  f  are  considered  to  be  electron  sources.  This 
description  is  possible  if  there  is  a  sharp  demarcation 
between  the  trapping  region  and  the  region  of  rapid 
atmospheric  scattering,  that  is,  the  bounce  loss  cone  is  well 
defined,  which  is  tme  forZ,  S;  1.2  [Walt,  1966], 

[9]  For  describing  the  stably  trapped  radiation  belt 
population  it  is  usual  to  introduce  bounce  and  drift 
averaging.  This  will  also  turn  out  to  be  useful  for  our 
numerical  solution  of  the  model  equation  (2).  The  bounce 
average  of  any  quantity  u  is  generally  defined  by  the  field 
line  integral 


(«)*=—/ — - — ds  (6) 

T A  /  V  cos  a 

taken  over  a  complete  bounce  cycle,  where  t b  is  the  bounce 
period  (defined  with  u  =  1).  The  bounce  averaged  transport 
equation  is  obtained  by  assuming  that  /  is  independent  of  s 
[. Retterer  et  al.,  1983]: 


where  t  is  time  and  Dao  is  the  diffusion  coefficient. 

[7]  It  is  useful  to  account  for  conservation  of  the  first 
adiabatic  invariant  by  a  change  of  variables  from  a  to  jc  = 
cos  oto,  where  oto  is  the  equatorial  pitch  angle.  With 
additional  terms  to  describe  the  radiation  belt  dynamics 
the  result  is 


df  df  df  cos  a  d  (  x 


where 


D  dl 

L'XX  o 

cos  a  ox 


+  Q- 


L, 

(2) 


cosa=  yfl  —  ^(1  -x2),  (3) 

v  is  the  constant  electron  speed,  5  is  the  distance  along  the 
magnetic  field  line,  wd  is  the  azimuthal  drift  rate,  4>  is  the 
azimuth  around  the  drift  shell  (or  drift  phase),  Dxx  =  (da J 
dx)2Daa  is  the  x  diffusion  coefficient,  Q  and  L  are  local 
source  and  loss  terms  respectively,  and  finally  B  and  B0  are 
respectively  the  local  and  equatorial  magnetic  field 
magnitude.  The  second  and  third  terms  account  respectively 
for  the  parallel  bounce  motion  [Davidson  and  Walt,  1977] 
and  the  azimuthal  drift  motion  [Selesnick  et  al.,  2003],  This 
transport  equation  neglects  radial  diffusion  and  energy 
diffusion,  which  are  significant  at  some  times  and  locations 
in  the  radiation  belt.  The  boundary  conditions  at  x  =  ±1  are 
dfldx  =  0  because  there  than  can  be  no  diffusion  into  or  out 
of  the  unphysical  x  region. 

[8]  The  source  and  loss  terms  are  included  in  (2)  to 
account  for  atmospheric  scattering  and  energy  degradation: 

Q  =  vcosa[5(i  -  4m]ffls+b(s  ~  4Jf£s\  (4) 

L  =  vcosa[6(s  —  s^tm)  +  b(s  —  (5) 


l+<“Al-is(w<0“>‘D  +  (e-I,‘  (7) 


which,  by  symmetry  about  x  -  0,  applies  to  the  interval  0  < 
x  <  1  with  dfldx  =  0  at  each  boundary.  The  bounce 
averaged  source  and  loss  terms  are  nonzero  only  at  x  values 
for  which  either  or  both  mirror  points  are  below  the  top  of 
the  atmosphere: 

(Q -L)b=±-  [($  -/)e(x"  - *)+(&  -/)©(*? -x)] ,  (8) 

'  '  b 


where  0  (x)  is  the  unit  step  function  (0  =  1  for  x  >  0,  0  =  0 
for  x  <  0).  The  x  values  xfl  and  correspond  to  electrons 
with  mirror  points  at  s  -  s^tm  and  ssatm  respectively.  These 
define  the  north  and  south  bounce  loss  cones. 

[10]  In  the  case  of  an  axisymmetric  magnetic  field,  such 
as  a  dipole,  it  is  simple  to  define  the  drift  average  as  the 
average  over  <j>.  Then  the  bounce  and  drift  averaged  trans¬ 
port  equation  is 

{Tbx{Dxx'>b’d~£)  +{Q~  L)b’d'  (9) 

where  (u)bid  is  the  bounce  and  drift  average  of  u.  This 
equation  does  not  strictly  apply  to  non-axisymmetric 
magnetic  fields  because  xb  and  iod  would  be  <j>  dependent, 
but  it  is  a  reasonable  approximation  for  the  Earth’s  inner 
magnetosphere. 

2.2.  Atmospheric  Backscattering 

[11]  The  north  and  south  backscattered  electron  distribu¬ 
tions  at  a  given  energy  E  are  related  to  the  distribution  / 
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entering  the  atmosphere  at  energies  greater  than  E,  because 
of  energy  loss  during  the  scattering  process: 

/"(a,£)=  fX  ['/(&,, o', E')S(a, E, a', £') sin a'do'dE' 

JE  JO 

(10) 

for  |  <  a  <  it  and 

pOO  pit 

fBss(a,E)  =  J  jf  /(^.a ',E')S(-k  -«,£,*  -  a 

•  sin  a'da'dE'  (11) 

for  0  <  a  <  where  a'  and  a  are  local  pitch  angles  at  s  = 
Sa,m  or  ssatm,  and  the  scattering  function  S  is  defined  for  a' 
values  entering  the  atmosphere  and  ot  values  leaving  the 
atmosphere  in  the  northern  hemisphere. 

[12]  The  scattering  function  S  is  the  backscattered  elec¬ 
tron  distribution  at  the  top  of  the  atmosphere  from  incident 
electrons  of  fixed  energy  and  pitch  angle  that  are  uniformly 
distributed  in  gyrophase.  To  calculate  S  we  use  the  EGSnrc 
Monte  Carlo  code  for  electron  and  photon  transport 
[Kawrakow  and  Rogers,  2000].  The  magnetic  field  is  that 
of  a  centered  dipole  at  L  =  3.8  and  200  km  altitude.  The 
incident  electrons  and  the  products  of  the  resulting  electron- 
7  ray  showers  are  followed  through  the  atmosphere  ac¬ 
counting  for  their  motion  in  the  local  magnetic  field. 
Electrons  of  each  sample  energy  and  pitch  angle  are  started 
at  an  altitude  of  200  km,  above  the  dense  atmosphere,  and 
the  emerging  electron  distribution  is  collected  at  205  km 
altitude  (slightly  higher  than  the  start  to  allow  for  some 
altitude  increase  during  the  initial  gyration  about  the  mag¬ 
netic  field).  The  atmosphere  below  200  km  is  constructed  of 
1  km  thick  layers  containing  atomic  and  molecular  O,  N, 
and  Ar  in  abundances  determined  by  the  NRLMSISE-00 
model  atmosphere  [Picone  et  al.,  2002]  averaged  over 
geographic  and  temporal  variations. 

[13]  Sample  results  from  this  calculation  are  shown  in 
Figure  1  for  two  different  incident  electron  energies  E' 
and  pitch  angles  a.  In  each  case  S  is  strongly  peaked  at 
energies  just  below  the  incident  energy,  but  there  is  also  a 
significant  contribution  to  the  total  number  of  backscat¬ 
tered  electrons  at  energies  E  well  below  the  incident 
energy  E1.  The  pitch  angles  of  the  backscattered  electrons 
are  generally  modified  significantly  from  those  that  would 
be  obtained  by  mirroring  of  the  incident  electrons  in  the 
dipole  field.  The  changes  in  backscattered  pitch  angle  and 
energy  increase  as  the  incident  pitch  angle  takes  the 
electrons  deeper  into  the  atmosphere.  The  total  albedo, 
or  fraction  of  incident  electrons  that  emerge  from  the 
atmosphere  at  any  energy  or  pitch  angle,  is  shown  in 
Figure  2  as  a  function  of  incident  pitch  angle  a'  for 
selected  values  of  the  incident  energy  E! .  The  total  albedo 
is  seen  to  decrease  significantly  with  increasing  energy 
except  at  the  higher  pitch  angles  where  the  reverse  is 
true.  This  is  because  only  the  lower-energy  electrons  are 
strongly  scattered  in  direction,  so  for  incident  trajectories 
with  mirror  points  well  below  the  top  of  the  atmosphere 
the  low  energies  can  easily  be  scattered  back  out  while 
the  high  energies  reach  deep  into  the  atmosphere  before 
losing  enough  energy  to  be  scattered.  For  incident  trajec¬ 


tories  with  mirror  points  in  the  high,  tenuous  atmosphere, 
the  low-energy  electrons  can  be  scattered  into  the  deep 
atmosphere  and  lost,  while  the  high-energy  electrons 
simply  mirror  without  scattering.  The  mirror  force  of 
the  magnetic  field  in  the  atmosphere  also  acts  to  increase 
the  albedo  relative  to  that  with  no  magnetic  field. 

2.3.  Diffusion  Coefficients 

[14]  Electron  pitch  angle  diffusion  by  cyclotron  resonance 
with  whistler  mode  waves  has  been  described  by  Lyons 
[1974]  and  Albert  [1999],  The  local  quasi-linear  diffusion 
coefficient  is  proportional  to  the  square  of  the  wave  mag¬ 
netic  field  bi,  and  is  a  complex  function  of  a  set  of  wave 
parameters  that  describe  the  distribution  of  wave  power  in 
frequency  and  propagation  direction,  and  of  the  local 
plasma  density.  At  L  =  3.5,  the  main  contribution  to  the 
whistler  mode  power  is  from  plasmaspheric  hiss.  Sample 
diffusion  coefficients  at  L  =  3.5,  calculated  by  the 
method  of  Albert  [1999],  are  shown  in  Figure  3  as  a 
function  of  x  for  two  different  electron  energies  and  for 
selected  magnetic  latitudes.  They  used  nominal  hiss  wave 
and  plasma  parameters  described  by  Albert  [1999],  based  on 
those  of  Abel  and  Thorne  [1998a],  For  each  magnetic 
latitude,  there  is  a  sharp  peak  at  low  x  that  is  a  result  of 
the  Landau  ( n  =  0)  resonance  of  electrons  near  their  mirror 
points  (a  near  90°).  The  continua  at  higher  x  values  are  due 
to  the  higher-order  resonances. 

[15]  The  calculated  model  diffusion  coefficients  can  be 
scaled  by  b2w  for  different  values  of  the  local  wave  power 
(Figure  3  used  bw  =  10  pT).  Separate  calculations  would  be 
required  if  any  of  the  other  wave  or  plasma  parameters  were 
varied,  but,  because  they  are  not  well  known,  we  use  just 
their  nominal  values.  The  sensitivity  of  the  diffusion  coef¬ 
ficients  to  the  various  wave  and  plasma  parameters  has  been 
described  by  Abel  and  Thorne  [1998b]  and  by  Albert 
[1999]. 

[16]  The  hiss  intensity  is  known  to  vary  with  magnetic 
local  time,  being  generally  highest  on  the  day  side  [Andre  et 
al.,  2002;  Selesnick  et  al.,  2003].  Therefore  we  have  chosen 
to  model  the  square  of  the  wave  field  b\.  as  varying 
sinusoidally  in  magnetic  local  time  with  a  maximum  at 
noon  and  a  minimum,  that  is  a  factor  of  10  lower,  at 
midnight.  The  intensity  is  also  known  to  increase  with 
geomagnetic  activity  [Smith  et  al.,  1974;  Andre  et  al., 
2002].  Therefore,  for  the  noon  maxima  we  have  chosen 
three  representative  values,  bw  =  31,  54,  and  100  pT,  to 
cover  a  range  of  magnetospheric  conditions.  We  assume  the 
wave  intensity  to  be  distributed  uniformly  along  the  mag¬ 
netic  field  lines. 

3.  Model  Solutions 
3.1.  Numerical  Method 

[17]  For  times  of  steady  radiation  belt  decay  the  time 
dependence  of  /  can  be  separated  from  the  spatial  depen¬ 
dence  as  an  exponential  with  a  time  constant  that  determines 
the  e-folding  electron  lifetime.  In  our  numerical  solution, 
the  finite  difference  approximation  to  the  transport  equation 
is  integrated  forward  in  time  from  an  initial  condition  until 
the  steady  decay  is  reached.  Then  the  lifetimes  can  be 
obtained  directly  from  the  decay  rate  of  the  numerical 
solution  or  from  an  s,  x,  and  4>  average  of  L  —  Q,  the  net 
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Figure  1.  Backscattered  electron  pitch  angle  distributions  at  205  km  altitude  from  electrons  entering  the 
atmosphere  with  energy  E!  and  pitch  angle  a!  for  three  energies  E  <  E',  calculated  with  the  EGSnrc 
simulation. 


loss  rate,  appropriately  weighted  by  the  spatial  dependence 

of  f 

[is]  The  distribution  function  /  adjusts  to  changing  con¬ 
ditions  on  three  timescales:  the  diffusion  timescale  1 

~  days,  the  drift  timescale  2rJ^d  ~  minutes,  and  the  bounce 
timescale  jb  ~  seconds.  A  numerical  solution  of  the  local 
transport  equation  (2)  over  all  three  timescales  is  impracti¬ 
cal,  because  the  time  step  must  be  a  small  fraction  of  the 
shortest  timescale.  Therefore  the  full  numerical  solution 
starts  with  the  bounce  and  drift  averaged  equation  (9)  to 
obtain  the  x  dependence  of /  on  the  diffusion  timescale  at  x 
values  outside  the  loss  cones,  proceeds  to  the  bounce 
averaged  equation  (7)  to  obtain  the  x  and  <j>  dependences 
of /in  and  near  the  drift  loss  cone  on  the  drift  timescale,  and 
finally  returns  to  the  local  equation  (2)  to  obtain  the  s,  x  and 
<t>  dependences  of/in  and  near  the  bounce  loss  cones  on  the 
bounce  timescale.  At  each  stage  the  solution  is  followed 
until  a  steady  exponential  decay  is  reached  after  several 
intervals  of  the  relevant  timescale  and  the  results  of  the  first 
two  stages  become  the  initial  condition  for  the  next  stage. 


[19]  The  pitch  angle  diffusion  operates  in  all  three  stages, 
but  in  the  second  and  third  stages  it  operates  over  smaller 
spatial  scales  determined  by  the  extent  of  the  drift  and 
bounce  loss  cones  that  are  filled  by  diffusion  during  a  single 
drift  and  bounce  period  respectively,  after  which  the  elec¬ 
trons  are  lost.  Therefore  the  drift  and  bounce  timescales  are 
sufficient  to  describe  the  diffusion  during  these  stages. 

[20]  The  main  computational  burden  occurs  in  the  final 
stage  where  the  full  solution  in  three  dimensions  ( s ,  <j>,  and 
x)  is  obtained.  The  initial  condition  for  this  solution  is 
improved  by  replacing  the  bounce  averaged  model  with  the 
backscattered  distribution  from  each  hemisphere  inside  the 
bounce  loss  cone  for  which  the  electrons  are  departing  that 
hemisphere.  This  introduces  the  asymmetry  about  x  =  0  that 
is  caused  by  the  differences  between  the  north  and  south 
bounce  loss  cone  angles  and  provides  an  initial  condition 
that  is  fairly  close  to  the  final  model. 

[21]  The  numerical  grid  in  s  and  x  at  a  single  grid  point  in 
<j>,  for  the  third  stage  of  the  numerical  solution,  is  illustrated 
in  Figure  4  (a  finer  grid  is  used  in  the  actual  calculation). 
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Figure  2.  Calculated  atmospheric  electron  total  albedo,  defined  as  the  total  number  of  backscattered 
electrons  at  any  energy  and  pitch  angle  as  a  fraction  of  the  number  of  incident  electrons  at  a  fixed  energy 
E '  and  pitch  angle  a,  versus  a'  for  selected  values  of  E! . 


The  grid  spacing  in  5  is  determined  by  converting  s  to  xm, 
the  x  value  for  electrons  mirroring  at  that  s  location,  and 
then  using  the  same  grid  spacing  for  xm  as  for  x  [Davidson 
and  Walt,  1977].  The  scales  in  Figure  4  are  nonlinear,  with 
significantly  finer  spacing  of  grid  points  near  the  atmo¬ 
spheric  boundaries  in  s  and  x  because  of  the  steeper 
gradients  there.  The  direction  of  the  bounce  motion  for  a 
given  x  value  determines  the  direction  of  the  advective 
integration  in  s.  When  the  integration  reaches  the  atmo¬ 
sphere  for  x  values  inside  either  bounce  loss  cone  then  the 
backscattered  distribution  is  calculated  from  (10)  or  (11) 
and  this  becomes  the  boundary  condition  for  integration  in  x 
over  the  next  half  of  the  bounce  motion.  For  x  values 
outside  both  bounce  loss  cones,  when  a  mirror  point  is 
reached  then  the  current  /  value  is  used  for  the  start  of  the 
next  half  of  the  bounce  motion. 

[22]  The  loss  cone  angles  are  calculated  from  the  IGRF 
magnetic  field  model.  However,  the  numerical  grid  along 
the  field  line  is  based  on  a  dipole  magnetic  field.  That  is,  the 
conversion  from  xm  back  to  s  is  based  on  the  dipole  formula 
for  B/B0,  the  ratio  of  local  to  equatorial  magnetic  field 
magnitudes.  This  is  a  good  approximation  for  the  inner 
magnetosphere,  particularly  given  the  weak  dependence  of 
the  model  solutions  on  s  (see  below).  A  dipole  approxima¬ 
tion  is  also  used  in  calculating  the  azimuthal  drift  rate  ud 
and  the  bounce  period  t*. 

[23]  The  backscattered  distribution  at  a  given  energy 
depends  on  the  /  values  at  higher  energies.  Therefore  it  is 
necessary  to  carry  out  the  numerical  solution  over  the  whole 


range  of  energies  simultaneously.  The  energy  spectrum  of 
the  stably  trapped  population  is  imposed  as  an  initial 
condition  for  each  stage  of  the  solution  by  a  normalization 
factor  at  x  =  0.  The  solution  is  sensitive  to  this  choice  of 
energy  spectrum  for  x  values  within  the  loss  cones. 

3.2.  Sample  Electron  Distributions 

[24]  Part  of  a  sample  numerical  solution  for  L  =  3.5  is 
shown  in  Figure  5.  It  includes  the  local  pitch  angle  a  and 
dipole  longitude  4>  dependences  of/at  a  fixed  energy  E  and 
at  a  location  5  along  the  field  lines  for  which  B/B0  =  40  (that 
would  correspond  to  an  altitude  of  1430  km  in  a  centered 
dipole  approximation  to  the  geomagnetic  field).  Figure  5 
illustrates  how  the  variation  with  ©  in  the  size  of  the  bounce 
loss  cones  causes  the  drift  loss  cone  intensity  to  vary  with 
longitude  as  a  result  of  the  balance  between  drift  and  pitch 
angle  diffusion  from  the  stable  trapping  region.  Then,  the 
resulting  intensity  in  the  drift  loss  cone  at  a  given  longitude 
determines  the  contents  of  the  bounce  loss  cones  at  that 
longitude.  The  longitude  <\>  is  measured  eastward  in  the 
dipole  magnetic  equator,  starting  from  0°  geographic  longi¬ 
tude,  and  increases  in  the  same  direction  as  the  electron  drift. 
Because  this  coordinate  system  rotates  with  the  Earth  and  the 
magnetic  local  time  dependence  of  the  diffusion  coefficient 
does  not,  it  is  necessary  to  perform  the  calculation  at  varying 
universal  times,  UT.  In  the  case  of  Figure  5,  UT  =  0  and  the 
noon  maximum  diffusion  coefficient  is  at  c|>  =  1 80°. 

[25]  Model  local  pitch  angle  distributions  at  two  fixed  <j> 
values  are  shown  for  selected  energies  in  Figure  6.  They 
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Figure  3.  Diffusion  coefficients  versus  x  =  cos  ao  for  energy  E  and  selected  magnetic  latitudes 
calculated  from  the  plasmaspheric  hiss  whistler  mode  wave  model  described  in  the  text.  See  color  version 
of  this  figure  in  the  HTML. 


give  a  more  detailed  view  of  the  bounce  loss  cone  distri¬ 
butions  over  a  range  of  energies.  Electrons  in  the  northern 
or  southern  bounce  loss  cone  entering  the  atmosphere  in  the 
corresponding  hemisphere  are  backscattered  to  form  the 
bounce  loss  cone  distributions  entering  the  opposite  hemi¬ 
sphere.  Pitch  angle  diffusion  during  the  bounce  motion 
increases  the  intensity  just  inside  the  bounce  loss  cones 
sufficiently  to  maintain  almost  the  same  backscattered 
distribution  on  the  next  bounce.  In  the  case  of  Figure  6a, 
for  example,  the  southern  bounce  loss  cone  is  larger  than  the 
northern  one  (because  of  the  south  Atlantic  anomaly,  or 
SAA)  so  that  the  electrons  diffusing  into  the  southern 
bounce  loss  cone  come  from  the  trapped  population,  where¬ 
as  those  diffusing  into  the  northern  bounce  loss  cone  come 
from  the  lower  southern  bounce  loss  cone  population. 
Therefore  the  backscattered  intensity  from  the  northern 
atmosphere  in  the  northern  bounce  loss  cone  (at  high  pitch 
angles)  is  much  lower  than  the  backscattered  intensity  from 
the  southern  atmosphere  in  the  northern  bounce  loss  cone 
(at  low  pitch  angles).  These  electrons  are  backscattered 


twice  per  bounce  period  whereas  those  outside  the  northern 
bounce  loss  cone  but  inside  the  southern  bounce  loss  cone 
are  backscattered  once  per  bounce  period,  from  the  southern 
atmosphere  only.  The  energy  spectra  in  the  bounce  loss 
cones  are  considerably  softer  than  those  in  the  trapped 
distribution.  This  is  a  result  of  the  generally  lower  albedo 
at  higher  energies  (Figure  2).  Also,  the  bounce  loss  cone 
spectra  in  Figure  6b  are  softer  than  those  of  Figure  6a 
because,  at  the  longitude  of  Figure  6b,  the  drift  loss  cone 
intensities  are  reduced  at  the  higher  energies  by  the  corre¬ 
spondingly  higher  drift  rates. 

[26]  The  sample  solutions  described  above  were  specifi¬ 
cally  at  the  northern  end  of  the  field  line.  The  solution  is 
different  at  the  southern  end  of  the  field  line  because  of  the 
diffusion  into  the  bounce  loss  cones  during  a  half  bounce 
period.  However,  in  these  examples,  the  gradient  djldx  is 
sufficiently  steep  near  the  bounce  loss  cone  boundaries  that 
the  differences  between  the  solutions  at  each  end  of  the  field 
line  would  not  be  apparent  on  the  scale  of  the  figures. 
Solutions  at  lower  energies  have  a  stronger  dependence  on  5 
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Figure  4.  Subset  of  the  numerical  simulation  grid  at  one  longitude.  The  resolution  is  coarser  than  in  the 
actual  calculations.  The  arrows  show  the  direction  of  the  electron  bounce  motion  for  either  side  of  oto  = 
90°  (pc  =  0). 
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because  of  the  longer  bounce  periods  [Davidson  and  Walt, 
1977], 

[27]  The  model  pitch  angle  distributions  well  inside  the 
bounce  loss  cones  at  low  and  high  pitch  angles  and  at  the 
higher  energies  show  significant  fluctuations  (Figure  6). 
These  are  caused  by  the  statistics  of  the  Monte  Carlo 
backscattering  calculation  (Figure  1)  and  do  not  adversely 
affect  the  final  model  simulations  because  they  only  occur 
where  the  backscattered  intensities  are  low. 

4.  Data  Comparisons 

[28]  We  have  selected  three  days  for  comparison  of  the 
model  results  with  data  taken  by  the  PET  instrument  [Cook 
et  al.,  1993;  Selesnick  et  al.,  2003]  on  the  low-altitude 
(~600  km),  polar  orbiting  satellite  SAMPEX.  These  days 
were  chosen  during  periods  of  relatively  steady  electron 
decay  following  magnetic  storm  injections.  They  were  also 
chosen  as  days  with  three  differing  average  geomagnetic 
activity  levels  and  for  which  there  is  little  evidence  of 
changes  in  the  pitch  angle  diffusion  rate  during  each  day.  The 
SAMPEX  orbit  crosses  L  =  3.5  ~60  times  per  day  at  varying 
longitudes  and  at  two  opposing  local  times.  The  data  set  from 
each  day  consists  of  the  electron  flux  at  each  crossing  from 
each  of  three  electron  counters  labeled  PI,  ELO,  and  EHI, 
that  measure  electrons  in  the  energy  ranges  of  0.6  MeV, 
~1.5  to  6  MeV,  and  ~2.5  to  14  MeV,  respectively.  Each 
data  point  is  at  a  given  set  of  B/B0,  (j>,  and  UT  values,  is  in 
the  northern  or  southern  hemisphere,  and  covers  a  distinct 


range  of  x  values  corresponding  to  electrons  with  mirror 
points  at  and  below  the  SAMPEX  altitude.  Equatorial 
pitch  angles  for  mirror  points  at  SAMPEX  vary  from 
~7°  to  ~9.5°.  Thus  a  particular  model  solution  is  used 
to  simulate  each  data  point  by  combining  the  model 
electron  distribution  function  at  each  satellite  location 
and  UT  of  the  measurement  with  the  detailed  instrumental 
energy  and  angular  response  functions  for  each  electron 
counter  (see  Selesnick  et  al.  [2003]  for  details  of  this 
method). 

[29]  The  models  for  each  day  differ  in  their  value  of  the 
noon  maximum  whistler  wave  field  bw  and  in  the  electron 
energy  spectrum  that  normalizes  the  electron  distribution  at 
x  =  0  ( ol0  =  90°).  All  other  model  parameters  are  the  same 
for  each  day.  The  normalizing  energy  spectra  are  from 
electron  data  taken  at  L  =  3.5  by  the  HIST  instrument 
[Blake  et  al.,  1995]  on  the  high-altitude  Polar  satellite  (see 
Selesnick  et  al.  [2003]  for  examples  of  these  spectra).  In 
reference  to  both  the  Polar  and  SAMPEX  data,  L  and  B 
values  were  calculated  with  the  IGRF  magnetic  field 
model. 

[30]  The  data  from  SAMPEX  and  the  model  simulations 
of  those  data  are  shown  in  Figure  7a  for  1997  day  186, 
Figure  7b  for  1998  day  151,  and  Figure  7c  for  1998  day 
159.  In  each  case,  the  observed  and  model  fluxes  for  each 
rate  counter  are  shown  as  a  function  of  dipole  longitude  cj>. 
Typically  the  highest  fluxes  are  observed  at  southern  lon¬ 
gitudes  from  ~300°  to  ~30°E  (the  SAA  region).  These  are 
stably  trapped  electrons,  with  mirror  points  that  do  not  reach 
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Figure  5.  Sample  simulated  low-altitude  electron  distribution,  showing  the  normalized  electron 
intensity  versus  local  pitch  angle  and  dipole  longitude  for  a  fixed  L  shell,  B/B0,  energy  E,  maximum  wave 
magnetic  field  bm  and  universal  time  (UT).  Dotted  (dashed)  curves  indicate  the  north  (south)  bounce  loss 
cone  angles.  See  color  version  of  this  figure  in  the  HTML. 


the  atmosphere  at  any  longitude.  The  lowest  fluxes  are  in  a 
similar  longitude  range  but  from  the  northern  hemisphere. 
These  are  untrapped  electrons  inside  the  bounce  loss  cone 
(BLC)  that  will  enter  the  atmosphere  within  the  next 
bounce.  Their  intensity  is  determined  by  the  pitch  angle 
diffusion  rate  into  the  bounce  loss  cones  and  by  the 
atmospheric  albedo.  In  the  longitude  range  from  ~30°  to 
~300°E,  the  fluxes  generally  are  at  intermediate  levels  and 
increase  toward  the  east.  These  are  the  quasi-trapped  elec¬ 
trons  in  the  drift  loss  cone  (DLC)  that  will  reach  the 
atmosphere  within  the  next  drift  orbit.  Their  intensity  is 
determined  by  the  balance  between  the  eastward  drift  and 


the  pitch  angle  diffusion  rate  from  the  stably  trapped 
population  into  the  drift  loss  cone. 

[31]  The  three  sample  days  show  similar  characteristics  in 
the  variation  of  the  electron  intensity  with  geographic 
location  of  the  data  points.  The  main  difference  between 
them  is  in  the  relative  intensity  levels  between  the  three 
regions,  those  of  the  stably  trapped,  quasi-trapped  (drift  loss 
cone)  and  untrapped  (bounce  loss  cone)  electrons.  In  the 
first  case,  1997  day  186,  the  pitch  angle  diffusion  rate  was 
lowest,  leading  to  less  filling  of  the  quasi-trapped  and 
untrapped  regions  relative  to  the  stable  trapped  region.  In 
the  third  case,  1998  day  159,  the  pitch  angle  diffusion  rate 
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Figure  6.  Simulated  low-altitude  pitch  angle  distributions  from  the  same  model  solution  as  in  Figure  5 
but  for  two  fixed  dipole  longitudes  ip  and  selected  electron  energies  E.  See  color  version  of  this  figure  in 
the  HTML. 


was  highest  and  so  the  quasi-trapped  and  untrapped  regions 
contain  relatively  high  electron  intensities.  The  second  case, 
1998  day  151,  was  intermediate  between  the  first  and  third. 
The  absolute  level  of  the  stably  trapped  intensity  also  varies 
somewhat  between  the  three  cases  as  a  result  of  the  prior 
history  of  the  radiation  belt.  This  determines  the  dynamic 
range  available  to  each  rate  counter  above  its  rate  of 
background  counts  (see  Selesnick  et  al.  [2003]  for  details). 

[32]  The  average  values  of  the  Dst  geomagnetic  index 
during  the  low,  medium,  and  high  diffusion  days  were  +1 
(1997  day  186),  -16  (1998  day  151),  and  -28  (1998  day 
159),  respectively,  so  they  were  all  relatively  quiet  days,  as 
is  typical  for  periods  of  radiation  belt  decay  following  a 
geomagnetic  storm.  The  correlation  between  the  Dst  index 


and  the  pitch  angle  diffusion  rate  (or  electron  lifetime)  is 
described  in  more  detail  by  Selesnick  et  al.  [2003]. 

5.  Discussion  and  Conclusion 

[33]  To  model  the  geographic  distribution  of  low-altitude 
radiation  belt  energetic  electrons  at  a  given  L  shell,  it  was 
necessary  to  include  both  the  bounce  and  drift  phase 
dependences  in  order  to  account  for  the  longitude  variations 
of  the  bounce  and  drift  loss  cone  electron  populations.  The 
model  included  pitch  angle  diffusion  at  a  rate  derived  from 
standard  models  of  whistler  mode  waves  in  the  plasma- 
spheric  hiss  frequency  band  and  of  the  plasma  density  at  L  = 
3.5.  Three  levels  of  wave  magnetic  field  intensity  were 


Figure  7.  Low-altitude  satellite  data  and  model  simulations  versus  dipole  longitude  from  the  three  electron  rate  counters 
(PI,  ELO,  and  EHI)  for  three  selected  days.  Each  data  point  and  its  corresponding  model  simulation  are  connected  by  a 
dashed  line.  Model  wave  fields  bw  correspond  to  (a)  low,  (b)  medium,  and  (c)  high  pitch  angle  diffusion  rates.  See  color 
version  of  this  figure  in  the  HTML. 
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chosen  to  normalize  the  diffusion  rates  for  comparison  of 
the  model  results  with  data  taken  on  days  with  three 
correspondingly  different  geomagnetic  activity  levels.  The 
model  also  included  electron  backscattering  from  a  model 
atmosphere  of  realistic  composition  calculated  with  the 
EGSnrc  Monte  Carlo  simulation. 

[34]  The  model  simulations  and  the  SAMPEX  electron 
data  are  in  reasonable  agreement  for  each  of  the  three 
sample  days  (Figure  7),  showing  that  the  model  includes 
the  main  factors  that  influence  the  low-altitude  electron 
distribution  during  times  of  relatively  steady  radiation  belt 
decay.  While  this  is  true  for  each  day  in  an  average  sense, 
there  are  instances  of  fairly  significant  discrepancy  between 
specific  data  points  and  their  simulated  values.  There  are 
several  possible  causes.  First,  the  model  parameters  were 
not  fit  to  the  data.  Only  the  wave  magnetic  field  bw  was 
varied  between  the  three  cases  and  this  was  simply  fixed 
at  the  three  values  judged  to  provide  examples  of  low, 
medium,  and  high  diffusion  rates  based  on  comparisons 
with  the  electron  data.  Second,  the  wave  parameters 
including  bw  and  those  affecting  the  wave  spectrum  and 
propagation  directions  are  not  well  known  even  in  an 
average  sense.  This  is  also  true  of  the  spatial  distribution 
of  the  waves  and  of  the  plasma  density  in  local  time  as  well 
as  along  the  magnetic  field  lines.  Third,  there  can  be 
changes  in  the  diffusion  rate  during  a  day  caused  by  varying 
wave  and  plasma  conditions.  While  the  sample  days  were 
chosen  to  minimize  this  possibility  by  comparison  with  data 
sets  obtained  on  adjacent  days,  there  must  always  be  some 
temporal  variability.  Fourth,  the  backscattering  calculation 
cannot  be  perfectly  accurate,  although,  even  with  the 
various  approximations  involved,  we  expect  that  it  should 
not  be  a  major  source  of  error  because  the  inputs  to  the 
calculation  are  relatively  well  known.  Fifth,  the  magnetic 
mapping  between  the  high-altitude  spectral  measurements 
on  Polar  that  are  used  to  normalize  the  model  and  the  low- 
altitude  locations  of  SAMPEX  may  introduce  some  error. 
Sixth,  the  model  neglects  pitch  angle  diffusion  by  lightning 
and  ground  based  VLF  generated  whistler  waves  and  other 
wave  modes.  While  these  are  significant  at  lower  L  values, 
their  impact  at  L  =  3.5  is  thought  to  be  minimal  in  a  bounce 
and  drift  averaged  sense  [Abel  and  Thorne,  1998a],  How¬ 
ever,  it  is  possible  that  the  lightning  and  VLF  signals  may 
make  some  significant  contributions  to  the  local  pitch  angle 
diffusion  coefficient  at  certain  times  and  longitudes  related 
to  the  occurrences  of  thunderstorms  and  the  geographic 
locations  of  VLF  transmitters.  These  could  be  added  to  the 
model  if  the  spatial  and  temporal  distribution  of  the  waves 
were  known.  Seventh  and  finally,  the  model  does  not 
include  radial  diffusion,  energy  diffusion,  and  possibly  other 
radiation  belt  processes  that  may  be  of  some  significance. 

[35]  While  there  is  a  long  list  of  possible  improvements,  it 
is  likely  that  the  greatest  source  of  error  is  in  the  plasma- 
spheric  hiss  wave  model  that  determines  the  pitch  angle 
diffusion  rate.  Some  evidence  of  this  can  be  seen  in  the 
comparisons  with  the  SAMPEX  data.  First,  the  higher 
energy  electron  measurements  (the  EHI  data)  tend  to 
disagree  with  the  simulations  to  a  greater  extent  than  those 
at  lower  energies.  This  may  be  indicative  of  an  inaccuracy 
in  the  energy  dependence  of  the  calculated  diffusion  coef¬ 


ficients.  Second,  the  normalization  of  the  model  by  the 
high-altitude  Polar  data  does  not  always  lead  to  an  accurate 
prediction  of  the  overall  stably  trapped  intensity  level 
observed  on  SAMPEX.  This  can  be  a  result  of  inaccurate 
distribution  of  wave  power  along  the  magnetic  field.  For 
example,  a  reduction  in  the  relative  wave  intensity  at  high 
magnetic  latitudes  would  reduce  the  equatorial  pitch  angle 
dependence  of  the  electron  intensity  in  the  stable  trapping 
region  and  increase  the  simulated  low-altitude  stably 
trapped  intensity.  These  considerations  show  that  improved 
measurements  of  the  spatial,  temporal,  and  frequency  dis¬ 
tribution  of  whistler  and  possibly  other  wave  modes  could 
allow  more  accurate  modeling  of  the  low-altitude  radiation 
belt  electron  distribution. 

[36]  Acknowledgments.  This  work  was  supported  by  NASA  under 
grant  NAG5-12766  and  by  the  Space  Vehicles  Directorate  of  the  Air  Force 
Research  Laboratory. 

[37]  Arthur  Richmond  thanks  Rumi  Nakamura  for  her  assistance  in 
evaluating  this  paper. 

References 

Abel,  B.,  and  R.  M.  Thome  (1998a),  Electron  scattering  loss  in  the  Earth’s 
inner  magnetosphere:  1.  Dominant  physical  processes,  J.  Geophys.  Res., 
103,  2385. 

Abel,  B.,  and  R.  M.  Thome  (1998b),  Electron  scattering  loss  in  the  Earth’s 
inner  magnetosphere:  2.  Sensitivity  to  model  parameters,  J.  Geophys. 
Res.,  103,  2397. 

Albert,  J.  M.  (1999),  Analysis  of  quasi-linear  diffusion  coefficients, 
J.  Geophys.  Res.,  104,  2429. 

Andre,  R.,  F.  Lefeuvre,  F.  Simonet,  and  U.  S.  Inan  (2002),  A  first  approach 
to  model  the  low  frequency  wave  activity  in  the  plasmasphere,  Ann. 
Geophys.,  20,  981. 

Blake,  J.  B.,  et  al.  (1995),  CEPPAD:  Comprehensive  energetic  particle  and 
pitch  angle  distribution  experiment  on  POLAR,  Space  Sci.  Rev.,  71,  531. 
Cook,  W.  R.,  et  al.  (1993),  PET:  A  proton/electron  telescope  for  studies  of 
magnetospheric,  solar,  and  galactic  particles,  IEEE  Trans.  Geosci.  Re¬ 
mote  Sens.,  31,  557. 

Davidson,  G.,  and  M.  Walt  (1977),  Loss  cone  distribution  of  radiation  belt 
electrons,  J.  Geophys.  Res.,  82,  48. 

Kawrakow,  I.,  and  D.  W.  O.  Rogers  (2000),  The  EGSnrc  code  system: 
Monte  Carlo  simulation  of  electron  and  photon  transport,  Tech.  Rep. 
PIRS-701,  Natl.  Res.  Counc.  of  Can.,  Ottawa,  available  at  http:// 
www.sao.nrc.ca/inms/irs/EGSnrc/pirs701/index.html. 

Lyons,  L.  R.  (1974),  Pitch  angle  and  energy  diffusion  coefficients  from 
resonant  interactions  with  ion-cyclotron  and  whistler  waves,  J.  Plasma 
Phys.,  12,  417. 

Picone,  J.  M.,  A.  E.  Hedin,  D.  P.  Drob,  and  A.  C.  Aikin  (2002), 
NRLMSISE-00  empirical  model  of  the  atmosphere:  Statistical  compari¬ 
sons  and  scientific  issues,  J.  Geophys.  Res.,  707(A12),  1468, 
doi:  1 0. 1 029/2002JA009430. 

Retterer,  J.  M.,  J.  R.  Jasperse,  and  T.  S.  Chang  (1983),  A  new  approach 
to  pitch  angle  scattering  in  the  magnetosphere,  J.  Geophys.  Res.,  88, 
201. 

Selesnick,  R.  S.,  J.  B.  Blake,  and  R.  A.  Mewaldt  (2003),  Atmospheric 
losses  of  radiation  belt  electrons,  J.  Geophys.  Res.,  108( A12),  1468, 
doi:  10.1 029/2003 JA01 0160. 

Smith,  E.  J.,  A.  M.  A.  Frandsen,  B.  T.  Tsurutani,  R.  M.  Thome,  and  K.  W. 
Chan  (1974),  Plasmaspheric  hiss  intensity  variations  during  magnetic 
storms,  J.  Geophys.  Res. ,  79,  2507. 

Walt,  M.  ( 1 966),  Loss  rates  of  trapped  electrons  by  atmospheric  collisions, 
in  Radiation  Trapped  in  the  Earth's  Magnetic  Field,  edited  by  B.  M. 
McCormac,  p.  37,  D.  Reidel,  Norwell,  Mass. 


J.  M.  Albert,  Air  Force  Research  Laboratory/VSBX,  29  Randolph  Rd., 
Hanscom  AFB,  MA  01731-3010,  USA. 

M.  D.  Looper,  The  Aerospace  Corporation,  P.O.  Box  92957-M2/260,  Los 
Angeles,  CA  90009-2957,  USA. 

R.  S.  Selesnick,  The  Aerospace  Corporation,  P.O.  Box  92957-M2/259, 
Los  Angeles,  CA  90009-2957,  USA.  (richard.s.selesnick@aero.org) 


11  of  11 


